Tuning the many-body interactions in a helical Luttinger liquid

In one-dimensional (1D) systems, electronic interactions lead to a breakdown of Fermi liquid theory and the formation of a Tomonaga-Luttinger Liquid (TLL). The strength of its many-body correlations can be quantified by a single dimensionless parameter, the Luttinger parameter K, characterising the competition between the electrons’ kinetic and electrostatic energies. Recently, signatures of a TLL have been reported for the topological edge states of quantum spin Hall (QSH) insulators, strictly 1D electronic structures with linear (Dirac) dispersion and spin-momentum locking. Here we show that the many-body interactions in such helical Luttinger Liquid can be effectively controlled by the edge state’s dielectric environment. This is reflected in a tunability of the Luttinger parameter K, distinct on different edges of the crystal, and extracted to high accuracy from the statistics of tunnelling spectra at tens of tunnelling points. The interplay of topology and many-body correlations in 1D helical systems has been suggested as a potential avenue towards realising non-Abelian parafermions.

Molecular-beam epitaxy (MBE) of monolayer 1T'-WTe 2 was carried out on highly oriented pyrolytic graphite (HOPG) and bilayer graphene (BLG) on 6H-SiC(0001) substrates. Prior to MBE, the HOPG substrates were mechanically cleaved in ambient conditions and then degassed at 300 ℃ overnight in ultrahigh vacuum (UHV) with a base pressure of 2 × 10 −10 mbar. For the epitaxial growth of BLG, 6H-SiC(0001) substrates were degassed at 600 ℃ overnight in UHV and then flash annealed at 2000 ℃ for a few cycles. Monolayer 1T'-WTe 2 was synthesised thereafter in an Omicron Lab10 MBE system (base pressure: 2 × 10 −10 mbar) with an electron beam evaporator (W source, purity = 99.998%) and a Knudsen diffusion cell (Te source, purity = 99.999%). We typically use a W:Te flux ratio of 1:280 and a substrate temperature of 165 ℃ for 2 hours to cover approximately 50% of the monolayer area.
A typical STM image of MBE-grown monolayer 1T'-WTe 2 island is shown in the left panel of Supplementary  Fig. 1a (i.e., the same island as shown in Fig. 1a of the main text). The edges of the island are generally disordered, with edges parallel or perpendicular to the atomic rows. For all our measurements, we chose regions with clean surface and well-defined edges as seen in the lower half of the island. The middle and right panels of Supplementary  Fig. 1a show the corresponding differential conductance dI/dV maps of the whole island, recorded at energies -70 meV and -30 meV. At -70 meV, both the bulk and edge local density of states (LDOS) Fig. 1: STM topography image and differential conductance dI/dV maps of a monolayer 1T'-WTe2. a Topography image of a monolayer 1T'-WTe2 island on HOPG (the same island in Fig. 1 of the main text) alongside its corresponding dI/dV maps at E − EF = −70 meV and −30 meV (recorded in constant height mode). The black squares correspond to the sample regions omitted from the dI/dV scanning to avoid the impurities depicted by the white spots in a. b, dI/dV map of the impurity-free lower half of the island in a, measured at energies from E − EF = −60 meV to E − EF = +50 meV. An enhanced LDOS persists along the sample edge for all energies, but the bulk is insulating between approximately -40 meV and +20 meV, with the LDOS matching that of the HOPG substrate. All the scale bars in this figure represent 5 nm.
substrate LDOS, indicating a metallic 2D bulk and edge. At -30 meV, inside the bulk energy gap (see also the insert to Supplementary Fig. 1f), only the LDOS at the edge is enhanced while the bulk LDOS is strongly suppressed. The presence of edge states surrounding an insulating bulk is consistent with the expected signatures of a quantum spin Hall (QSH) insulator. Further dI/dV maps at different energies are shown in Supplementary Fig. 1b. We find that the enhanced edge state LDOS persist at all energies, while measurable bulk LDOS is only observed between ≃ −40 meV and ≃ 30 meV.

S2. TEMPERATURE-DEPENDENT SCANNING TUNNELLING SPECTROSCOPY OF THE TLL
The tunnelling conductance dI dV (E, T ) at finite temperatures is modelled by the convolution of the TLL tunnelling density of states ρ TLL and the first derivative of the Fermi-Dirac distribution f (E), i.e., where   Supplementary Fig. 2 shows a full set of data, demonstrating universal scaling for different edge terminations on both BLG and HOPG. For each substrate and edge, the tunnelling spectra collapse onto a single universal curve that fits to Supplementary Equation (2) with a common α. We find α = 0.45 ± 0.15 (K = 0.40 ± 0.05) for WTe 2 /HOPG Y-edge, α = 1.05 ± 0.15 (K = 0.26 ± 0.02) for WTe 2 /HOPG X-edge, α = 0.98 ± 0.15 (K = 0.27 ± 0.03) for WTe 2 /BLG Y-edge, and α = 1.35 ± 0.15 (K = 0.22 ± 0.01) for WTe 2 /BLG X-edge. All values extracted for K fall into the normal distributions shown in Fig. 3 of the main text.

S3. THE DEPENDENCE OF THE ZBA CONDUCTANCE ON THE LOCK-IN AMPLITUDE
All scanning tunnelling spectroscopy was carried out using standard lock-in techniques with a modulation amplitude of V lock−in = 2 mV. This modulation amplitude has been optimised to maximize the signal-to-noise ratio, while minimizing Gaussian broadening of the instrumental response. Supplementary Fig. 3a shows our amplitude calibration, comparing spectra of the TLL pseudogap, measured at T = 4.5 K and with lock-in amplitudes ranging from V lock−in = 0.25 mV to 10 mV. Spectra taken at V lock−in < 3 mV fit well to the TLL model with a common α = 0.83 (K = 0.3). Instrumental broadening has a measurable effect only at V lock−in > 3 mV, more clearly seen in the zero-bias (E = E F ) tunneling conductance ( Supplementary Fig. 3b). For spectra taken at V lock−in < 1 mV, we find that the signal-to-noise ratio substantially degrades. An AC modulation of 2 mV therefore provides a balance between optimized signal-to-noise ratio and suppressed instrumental broadening. b, Fitting residue calculated for the traces in a to estimate the χ 2 value. We find acceptable fits to a CG model only for Ω = (4 ± 2) meV. None of the CG fits can model the data over the full energy range, but selectively only at the ZBA minimum or its shoulder.

S4. RULING OUT A COULOMB GAP AS ORIGIN OF THE ZERO BIAS ANOMALY
An alternative mechanism, sometimes invoked to explain the presence of a zero bias anomaly (ZBA) is the formation of an Efros-Shklovskii Coulomb gap (CG) [1,2], where long-range interactions between electrons cause the density of states (DOS) near the Fermi level to be suppressed. To rule out the presence of a CG, we fit the data in Fig. 2 of main text to the CG model. For the quasi-1D case, the DOS near the Fermi level is expressed as [3,4]:  where E, E ′ represent energy, is the reduced Planck constant, t is time, ρ 0 is a proportionality factor, and Ω is the frequency scale which is defined by: Here, f 0 represents the effective electron-electron interaction potential corresponding to the long-range interactions between electrons. At T → 0 K, f 0 is equal to the interaction potential W (q = 0), only when the exchange part of the interactions is considered. D 0 is the diffusion coefficient. This makes Supplementary Equation (3) a single-parameter fit at each temperature with the frequency scale Ω the only variable.
Using the data presented in Fig. 1g of the main text, we illustrate how we differentiate between TLL and CG models. The same data, reproduced in Supplementary Fig. 4a, was spatially averaged to achieved improved noise levels and fit to Supplementary Equation (2) and (3), respectively. The corresponding fitting residues are shown in Supplementary Fig. 4b. We find that acceptable fits to a CG model can only be obtained for Ω values ranging between 2-6 meV. However, the fit can only reproduce, either, the minimum of the ZBA or the shape of its shoulder, but never both at the same time. In stark contrast the data well-described over the full energy range within the QSH gap by a TLL model with α = 0.575 ± 0.025. This is further reflected in a lower least squares value χ 2 = 1.75217 for the TLL fit, compared to the best CG fit (χ 2 = 5.16056 at Ω = 3 meV for this data set). The same procedure has subsequently been applied to all universal scaling data sets on both HOPG and BLG substrates, including Supplementary Fig. 5a,b and Fig. 2e,f of the main text to determine the respective Ω at T = 4.5 K. The temperature dependence of the fits to both TLL (black lines) and CG models (magenta lines) between 4.5 K ≤ T ≤ 25 K are shown in Supplementary  Fig. 5c,d. We find that the TLL model (black lines) provides better agreement with our experimental data for all temperatures and both substrates, with the differences between CG and TLL most clearly visible in the semi-log plot of the zero-bias conductance shown in Fig. 2e,f of the main text. Correspondingly, Supplementary Fig. 5e shows the temperature dependence of the QSH bulk gap, which persists up to T ≈ 25 K.

S5. STATISTICAL DISTRIBUTIONS FOR K AND α
The histograms in Supplementary Fig. 3a,b in the main text show the Luttinger parameter K extracted at different edges of the 1T'-WTe 2 crystal on HOPG and BLG substrates, respectively. A total of 69 locations were probed across different samples on both substrates and with 5 different tunneling tips. All four histograms fit well to Gaussian normal distributions, from which we extract mean and standard deviations for the Luttinger parameters K as listed in Table 1. The extracted values have subsequently been plotted as a function of substrate relative permittivity in Fig. 3c of the main text.
To further confirm normality of the four distributions, we produce Q-Q plots for all data sets, shown in Supplementary Fig. 6. For perfect normal distributions, the theoretically expected quantiles (which are obtained by taking the inverse of the cumulative distribution function of the Gaussian distribution) versus the raw data quantiles will lie on the identity line (red solid line) as shown. Slight skew (deviation of data quantiles from the identity line) is only observed for the X-edges of WTe 2 /HOPG and WTe 2 /BLG. Otherwise, the data is well described by normal distributions.  Fig. 6: Q-Q plots for the histograms in Fig. 3a,b in the main text. We find all four data sets, when plotted against the expected values, lie on the identity line, confirming normality of the distributions. Slight deviations (skew) exists only for the X-edges of WTe2/HOPG and WTe2/BLG.
Given that the distributions are overlapping and that the sample number is relatively small, we further confirm that all distribution are statistically distinct. To this end, we run two-sample T −tests, applicable for small populations, with the null hypothesis that any pair of the four distributions are part of the same wider distribution. This procedure yields a p-value quantifying the statistical significance of the difference in mean.  Fig. 3a,b of the main text by fits to Gaussian normal distributions. The mean values of K, standard deviation σ, and standard error of the mean SEσ = σ √ n , where n is the number of data points, are listed here and incorporated in Fig. 3c Fig. 3a,b of the main text. The test returns p = 1 if any pair of distributions are identical, and p = 0 if they are perfectly distinct. For all pairs of distributions, we find p < 0.05, implying that the distributions are statistically distinct within a confidence interval of 95%.

Type
Type BLG X-edge BLG Y-edge HOPG X-edge HOPG Y-edge BLG The T −value is defined as, where µ 1,2 and σ 1,2 are the sample mean and standard deviation of distributions 1 and 2, respectively, while m and n are their number of samples. The degree of freedom df of the T −test (generally for unequal variances) represents the effective number of independent variables in the two distributions, defined as Once the T −value and df are obtained, the p−value can be calculated as where cdf is the cumulative distribution function for the Student's T −distribution. Our T −test yields p−values p < 0.05 for any pair of two distributions (see Table 2), confirming that that each distribution is statistically distinct with a statistical significance of > 95%.

S6. ELECTROSTATIC MODEL TO ESTIMATE THE LUTTINGER PARAMETER K
For a theoretical estimate of the Luttinger parameter K, we calculate the electrostatic energy stored inside a threedimensional space-charge distribution along the sample edge. Such model was first proposed by Teo et al [5] and recently successfully employed by Stühler et al [6]. The charge density can be expressed as, where L is the edge length, w = 0.7 nm is the vertical thickness of the charge distribution (which is given by the 1T'-WTe 2 monolayer thickness), and Q is the total charge. The charge distribution is concentrated along the sample edge-vacuum interface, and decays exponentially into the 2D bulk with a characteristic decay length ξ ≃ 2 nm (see Fig. 1 of the main text). To model dielectric screening by a substrate, we consider a substrate-dependent mirror charge density at a vertical separation of 0.5 nm, corresponding to the van-der-Waals gap between the 1T'-WTe 2 monolayer and its substrate [7]. The image charge distribution is given by: where ǫ r is the dielectric constant of the substrate. The electrostatic energy stored in the edge state-image charge system is then evaluated numerically, via: where r, r ′ are coordinates of the individual point charges. The first term in Supplementary Equation ( 10) represents the unscreened interaction of real charges at the edge, while the second term represents the electrostatic interaction with the image charge. For an isotropic Coulomb interaction in the long-wavelength limit (q = 0), the Coulomb interaction strength W is equal to the potential W = W (q = 0) = U E L [5,6], which we then use to obtain K (Eq. 1 of the main text). We note that our approach differs from that of Stühler et al [6], in that we evaluate Supplementary Equation( 10) numerically rather than approximating it by a logarithmic interpolation between the two limits ξ ≪ w ≪ L (i.e., negligible decay length) and w ≪ ξ ≪ L (i.e., negligible channel width). Similar to bismuthene, the edge state decay length in 1T'-WTe 2 is comparable to the thickness of the atomic monolayer (ξ ∼ w) [7,8]. Hence, neither limit strictly applies to our samples. We find that a logarithmic interpolation of W , where: provides reasonable agreement only at small ǫ r (as in the case of Ref. [6]). However, this interpolation fails in the perfect metal limit (ǫ r → ∞), where it predicts that K → 1, which corresponds to the case of non-interacting electrons. Indeed, the finite interaction strength observed in our samples (K ∼ 0.35), even in the metallic limit, is more accurately represented by a full numerical evaluation of Supplementary Equation ( 10), with the van-der-Waals gap between the sample and substrate taken into account. The agreement of this space-charge model with the experimental data is further confirmed by tight-binding calculations (see Fig. 3c of the main text and Sec. S7).

S7. FULL TIGHT-BINDING MODEL TO CALCULATE THE LUTTINGER PARAMETER K
The electronic structure in 1T'-WTe 2 was modeled by constructing a real space eight-band Hamiltonian for two Te and two W atoms in a rectangular unit cell. This low-energy Hamiltonian contains contribution from Te p x orbitals and W d x 2 −y 2 orbitals, which reflects the distorted tetrahedral crystal structure of the monolayer [9][10][11]. Band structures from this model are shown in Fig. 3d,e of the main text for the Y-and X-directions of the Brillouin zone. In addition, we include the hopping matrix elements with the Rashba spin-orbit coupling terms. An example of a real-space tight-binding calculation of the LDOS for realistic edge disorder is shown in Supplementary Fig. 7, showing an STM topography image (a) of the same X-edge as in Fig. 1 and Supplementary Fig. 1, as well as with extracted roughness overlaid (e). Supplementary Fig. 7(b-d) and (f-h) compare measured dI/dV maps with calculated LDOS maps at energies indicated. The LDOS has been calculated for the particular atomic configuration at the edge as extracted / indicated in (e), and shows good agreement with the modulations observed in the experiments.
To obtain estimates for the Luttinger parameter, we construct a mean-field tight-binding Hamiltonian for 1T'-WTe 2 , which has the following form: where t is the hopping integral, indices i, j refer to the unit cell, µ, ν = (1, 2, 3, 4) represent the atoms or orbitals within each unit cell, and h.c. stands for hermitian conjugate. The operators c † iµσ and c jνσ ′ are the creation and annihilation operators at a site with orbital and spin quantum numbers σ = (↑,↓). Using the Bogoliubov-de Gennes (BdG) transformation, these operators can be written as: x y a 30meV 10meV  Supplementary Equation (14)), we self-consistently calculate the electron density at each lattice site. The mean-field charge distribution ρ iµσ is given by: where f (E n ) is the Fermi-Dirac distribution at energy E n . The mirror charge density ρ m is then calculated using Eqs. 15 and 9 for 0 ≤ ǫ r ≤ 1000, with a 0.5 nm separation between the substrate and WTe 2 . By evaluating the U E (Supplementary Equation (10)) using the charge densities from tight-binding calculations for different values of ǫ r , we obtain K using Eq. 1 of the main text. An important aspect to consider for the validity of our formalism is the assumption of g 2 = g 4 in the perturbative limit, where g 2,4 are density-density interaction terms, corresponding to the Hamiltonian contributions, where ψ R↑ annihilates a right mover electron with spin up, and ψ L↓ annihilates a left mover electron with spin down. We note that the assumption g 2 = g 4 is reasonable if the range of the Coulomb interaction is significantly extended. As shown in Supplementary Fig. 8, the Coulomb interaction has significant contributions up to the fifth nearest neighbour, owing to the finite width of the helical edge channel (note that the non-monotonous behaviour is due to the variation in electron density between the Te and W ions). We anticipate that our extraction of the Coulomb interaction strength may stimulate further in-depth theoretical treatment, for instance using density-matrix renormalization group (DMRG) methods. Lattice position Coulomb Interaction (eV) 6 8 10 Supplementary Fig. 8: Calculated Coulomb interaction as a function of the lattice position. The Coulomb interaction has significant contributions up to the fifth nearest neighbour, corresponding to the fact that the helical edge channel has a finite width.